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1 Introduction 

Factorization of ordinary and partial linear differential operators (LODOs and LPDOs) is a 
very well-studied problem and a lot of pure existence theorems are known. For LODOs it is 
proven that a factorization is unique while for LPDOs even uniqueness is not true any more 
and in fact parametric families of factorizations can be constructed for a given LPDO as will be 
demonstrated below. First constructive method of factoring second order hyperbolic LPDO 
in the form 

d^dy + adx + hdy + c, (1) 

belong to Laplace who formulated it in terms of invariants a = c — ah — ax and b = c — ah — hy 
now called Laplace invariants. An operator is factorizable if at least one of its Laplace in- 
variants is equal to zero. Various algorithms are known now for factoring LPDOs over different 
differential fields beginning with the simplest field of rational functions [T] . 

Recently two papers (|2], 0) on factoring arbitrary order LPDOs have been published. In 
[2j a modification of well-known Hensel lifting algorithm (see, for instance, jlj) is presented 
and sufficient conditions for the existence of intersection of principal ideals are given. These 
results are applied then to re-formulate the factorization formulae for second and third order 
operators from the ring D = Q{x,y)[dx,dy] obtained by Miller (1932) in terms of principal 
intersections. Authors do not claim to construct general algorithm of factoring arbitrary order 
LPDOs but say that this approach is applied to "execute a complete analysis of factoring and 
solving a second order operator in two variables. Some results on factoring third-order operators 
are exhibited" . On the other hand, this approach can probably be generalized to many variables. 



In jS] necessary and sufficient conditions are given for factoring of bivariate LPDOs 
of arbitrary order with coefficients being arbitrary smooth functions. In jH] it was shown that 
this procedure called now BK-factorization principally can not be generalized on the case of 
more than two variables. In was also shown that conditions of factorization found in [H] are 
invariants under gauge transformation and classical Laplace invariants are particular case of 
this generalized invariants. In this paper we re-formulate BK-factorization in more suitable for 
applications invariant form and illustrate it with a few examples, give a sample of the symbolical 
implementation of this method in MATHEMATICA and also discuss some possibilities to use 
this method for approximate factorization of LPDOs. 
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2 Helsel descent and BK-factorization 



We begin this paper with a brief comparison of Hensel descent and BK-factorizations in order 
to show merits and draw-backs of each method, just for comphcity of presentation. 

• In 13] factorization of a bivariate LPDO is looked for in the form 

An= ""^^^A = (P^^^ + P^^y +P^){Y1 P^k^'^^l) (2) 

j+k<n j+k<n 

and in P] - operator of m > 2 independent variables is regarded and for m = 2 factoriza- 
tion is looked for in the form 

An= Yl ^^^^i^'y = ( E P^^^'A) ( E P^'^^'^i) ^'^^ l + r = n. (3) 

j+k<n j+k<l j+k<r 

• In jS] coefficients ajk are arbitrary smooth functions, for instance trigonometric functions; 
in [2] conditions for reducibility of an operator are studied when "coefficients are from 
a universal field of zero characteristic" , while "studying factorization algorithms we 
will assume that the input operators are from the ring Q{xi, Xm)[di, dm]" This 
suggestion is necessary: "From now on the coefficients of a given second-order operator 
are assumed to be from the base field Q{x,y). This is necessary if the goal is to obtain 
constructive answers allowing to factorize large classes of operators" ([2], Sec. 3); "In this 
section we study third-order LPDOs from the ring D = y)[dx, dy]." (|^, Sec. 4). 

• In |3] it was shown that in generic case factorization can be constructed explicitly and 
algebraically, while in (Sec. 5) it is concluded that "the factorization problem for second- 
and third-order differential operators in two variables has been shown to require the 
solution of a partial Riccati equation, which in turn requires to solve a general first-order 
ODE and possibly ordinary Riccati equation. The bottleneck for designing a factorization 
algorithm for a LPDO is general first-order ODE which make the full problem intractable 
at present because in general there are no solution algorithm available." 

• In |2] it is pointed out that BK-factorization procedure has to be modified in some way 
(presently unknown to authors) in order to proceed with operators of more than 2 inde- 
pendent variables, while in [2j (Sec. 5) it is written that "some of the results described in 
this article may be generalized to any number of independent variables." 

• In IHj it shown that BK-factorization procedure gives rise to construction of the whole class 
of generalized invariants particular case of them being classical Laplace invariants. This 
leads to a possibility to factorize simultaneously the whole class of operators equivalent 
under gauge transformation (see next Section) while Helsel descent is used for factoring 
of a one specific operator. 

We summarize all this in the Table below. 



Property / Method 


BK-factorization 


Hensel descent 


Order of operator 


n 


n 


Operator's coefficients 


arbitrary smooth functions 


rational functions 


Number of variables 


2 


possibly > 2 


Conditions of factorization 


necessary and sufficient 


sufficient 


Form of factorization 


factors of order 1 and {n — 1) 


factors of order k and {n — k) 


Formulation in terms of 


explicit formulae for fact, coefficients 


existence of ideals intersection 
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In the next Sections we demonstrate some other interesting properties of BK-factorization - 
first of all, that it has invariant form and can be used therefore to factorize simultaneously the 
whole classes of equivalent LPDOs. Second, the use of this invariant form of BK-factorization 
for construction of approximate factorization for LPDEs to be solved numerically. 



3 Invariant Formulation 

We present here briefly main ideas presented in j3j , IS] beginning with the definition of equivalent 
operators. 

Definition The operators A, A are called equivalent if there is a gauge transformation that 
takes one to the other: 

Ag = e~'^Aie-^g). 

BK-factorization is then pure algebraic procedure which allows to to construct explicitly a 
factorization of an arbitrary order LPDO A in the form 

A:= J2 ajk&y = Lo V.k^A 

j+k<n j+k<{n-l) 

with first-order operator L = dx — 0dy + p where is an arbitrary simple root of the 
characteristic polynomial 

n 

V{t)=J2an-k,kr-\ V{0) = O. (4) 

k=0 

Factorization is possible then for each simple root of (^J iff 
for n = 2 ^ ^2 = 0, 
for n = 3 ^ = 0, & /31 = 0, 
for n = 4 ^ /4 = 0, & /41 = 0, & U2 = 0, 

and so on. All functions I2, h, hi, U Ui, U2, ■■■ are explicit functions of aij and 0. 



Theorem All /2, /s, /31, .... are invariants under gauge transformations. 



Definition Invariants I2, 13, /31, .... are called generalized invariants of a bivariate operator of 
arbitrary order. 

In particular case of the operator its generalized invariants coincide with Laplace invari- 
ants. 



Corollary If an operator A is factorizable, then all operators equivalent to it, are also fac- 
torizable. 



As the first step of BK-factorization, coefficients pij are computed as solutions of some 
system of algebraic equations. At the second step, equality to zero of all generalized invariants 
lij = has to be checked so that no differential equations are to be solved in generic case. 
Generic case corresponds to a simple root of characteristic polynomial, and each simple root 
generates corresponding factorization. Moreover, putting some restrictions on the coefficients 
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of tlie initial LPDO a^j as functions of x and one can describe all factorizable operators 
in a given class of functions (see Example 5.3 in jSj). The same keeps true for all operators 
equivalent to a given one. Equivalent operators are easy to compute: 

e"^9^e'^ = dx + Lfx, e^^dye"^ = dy + Lpy, 

e-'^d^dye^ = e-^d.e^e-'^dye^ = {d, + ip,) o {dy + ipy) 

and so on. Some examples: 

• = d,dy + xd, + l = d,{dy + x), k{A,) = 1-1-0 = 0; 

• A2 = d^dy + xd^ + dy + x + l, A2 = e'^'Aie''; l2{A2) = {x + 1) - 1 - x = 0; 

• A3 = d,dy + 2xd^ + {y + l)dy + 2{xy + x + l), A3 = e-^Mae"?'; kiA^) = 2{x + 1 + 
xy) - 2 - 2x{y + 1) = 0; 

• Ai = d^dy + xd^ + (cosx + l)dy + xcosx + x + 1, A^ = e-^'^^Ase"'"^; hiA^) = 0. 

Generic case which can be treated pure algebraically by BK-factorization corresponds to a 
simple root of characteristic polynomial. Each multiple root leads to necessity of solving 
some Ricatti equation(s) (RE). If appeared RE happens to be solvable, such a root generates 
a parametric family of factorizations for a given operator. For instance, well-known Landau 
operator 

C + ^dl^ + + (2x + 2)dl + d, + i2 + x)dy 

has characteristic polynomial with one distinct root 0i = —x and one double root 02,3 = 0. 
Factorization then has form 

{d^ + r){d^-r + 2){d^ + xdy) 
where r is a solution of Ricatti equation 

l-2r + a^.(r) + r2 = 

which is easily solvable: 

1 

r = 1 + 



x + Y{y) 

with arbitrary smooth function Y{y) of one variable y so that factorization has form 

A={d^ + l+ , \. J (9. + 1 + ^^s/)- 

x + Yiy) x + Y{y) 

Notice that to factorize an ordinary differential operator it is always necessary to solve 
some RE. Nevertheless, just formal application of BK-factorization will produce all the linear 
factors in the case when corresponding RE are solvable. For instance, the factorization has 
been constructed in [Oj 

2 x'^ — 1 y/2 — 1 
xd^^^ + (x^ - l)d^^ - xd^ + — -1 = {d^-\ ){xd^ - V2)((9^ H ). 

while both RE appearing at the intermediate steps are solvable. 

These two last examples show the main difference between factorizing of ordinary and 
partial differential operators - LODO has always unique factorization while LPDO may has 
many. An interesting question here would be to compute the exact number of all possible 
factorizations of a given LPDO into all linear factors (its upper bound is, of course, trivial: 
n\). A really challenging task in this context would be to describe some additional conditions 
on the coefficients of an initial operator which lead to solvable RE. 
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4 Left and Right Factors 



Factorization of an operator is the first step on the way of solving corresponding equation. But 
for solution we need right factors and BK-factorization constructs left factors which are easy 
to construct. On the other hand, the existence of a certain right factor of a LPDO is equivalent 
to the existence of a corresponding left factor of the transpose of that operator. Moreover 
taking transposes is trivial algebraically, so there is also nothing lost from the point of view of 
algorithmic computation. 



Definition The transpose A* of an operator A = ^ a^d"", 5" = ■ ■ ■ 5°". is defined as 

A*« = ^(_l)Ha"(a„w). 

and the identity 



implies that 



Now the coefficients are 



with a standard convention for binomial coefficients in several variables, e.g. in two variables 

a\ _ /(ai,a2)\ _ /aA /'a2 

In particular, for order 2 in two variables the coefficients are 

ttjk = ttjk, j + k = 2; flio = — flio + '2dri.a2o + dyUu, agi = — Oqi + d^an + 2dyao2, 

o-oo = floo ~ dxttio — dyttoi + dla2o + dxdxttu + dyao2- 
For instance, the operator 

dxx - dyy + ydx + xdy + ^{y^ - x^) - 1 (5) 

is factorizable as 

[dx + dy + liv-x]] [...] 
and its transpose A\ is factorizable then as 

[...] [dx-dy+'^{y + x)]. 

Implementation of the BK-factorization for bivariate operators of order n < 4 is therefore 
quite straightforward and has been done in MATHEMATICA while all roots of characteristic 
polynomial are known in radicals. For instance, for the operator (0) with 2 simple roots we get 
one factorization 

[dx -dy + \{y + x)] [dx + dy + \{y - x) 
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corresponding to the first root while in the case of the second root, generalized invariant is 
equal to 2. 

If n > 5 the problem is generally not solvable in radicals and very simple example of non- 
solvable case is: — 4x — 2 = 0. Thus, to find solutions in radical for n > 4 one needs some 
constructive procedure of finding solvable Galois group but this lies beyond the scope of the 
present paper. 

5 Approximate Factorization 

An interesting possible application of the invariant form of BK-factorization is to use it for 
construction of approximate factorization of a given LPDE, in the case when exact factorization 
of corresponding LPDO does not exists. Indeed, as a results of BK-factorization one gets 

(1) factorization coefficients {p'fj} for the z-th factorization of a given operator, and 

(2) generalized inv£triants with all pfj, being explicit functions of the coeffi- 
cients of initial operator Ojj . 

In numerical simulations coefficients a.y of the equation are always given with some non-zero 
accuracy, say e > 0, which means that it is enough to construct an approximate factorization in 
the following sense. One has to find restrictions on the coefficients aij of an initial LPDO which 
provide < e with a given accuracy < £ << 1. Many different strategies are possible 
here, we just give a brief sketch of two approaches we are working on right now: 

5.1 Quantifier Elimination 

We illustrate this idea on the simple example of a hyperbolic operator d^x — dyy + cii^d^ + aQidy + 
aoo with hnear polynomial coefficients. What we have is: 

aoo(x, y) = hx + b2y + h, 
aw{x,y) = csx + C2y + ci, 
aoi(x, y) = dsx + day + di, 
a function constructed from general invariants 

^ _ ^3 - S2 {SsX + S2y + Sif 

2 ^ 4 

with Si = Ci — di. 

What we need is: 

To find some function(s) F — F{aij) such that if F{aij) = 0, then 

— £ < aoo — < for some constant < e << 1, 

i.e. to find some conditions on the initial polynomials which provide that function TZ differs 
not too much from one these polynomials, namely aoo- 

Notice that simple symmetry considerations allowed us to reduce number of variables needed 
for CAD calculations. Initially we had 9 variables 63, 62, 6i, C3, C2, Ci, (is, c^i but in fact it is 
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enough to regard only 6 variables 51,52,53,61,62,63. Nevertheless, the computation time may 
become crucial while using this approach due to the substantial number of variables. On the 
other hand, this approach allows us work generally on the operator level including initial and/or 
boundary conditions first at some later stage. 



5.2 Auxiliary Operator 

Another approach is to construct a new auxiliary operator with coefficients ciij = f{x,y))aij 
for all or for some of the coefficients Uij of the initial operator, keep invariants (almost) equal 
to zero and find function(s) f{x, y) minimizing the differences between the coefficients of initial 
and new operators. In this way an auxiliary operator is constructed which can be regarded as 
an approximate operator for the initial operator. Of course, it does not mean that solutions of 
the initial and approximate operators will be also close but simple properties of linear opera- 
tors show that it is necessary (but not sufficient!) step on the way of construction of a good 
approximate solution of a given LPDE - in the case of a well-posed problem, of course. In 
particularly, it means that one have to introduce proper metrics in the space of operators and 
in the space of solutions. Choice of the both metrics and of a function / will depend on (1) 
coefficients of the initial operator; (2) class of functions in which we are looking for a solution; 
(3) initial and/or boundary conditions. 

To demonstrate all this let us regard two different un-factorizable modifications of the 
operator (0): 

A = dxx- dyy + ydx + xdy + -(2/^ - x^)-l (6) 
with kiA) = liy"^ - x^) and 

B = dxx — dyy + sin y dx + cos x dy + 2 (^^^^ y ~ '^^^^ (7) 

with hiB) = I (cos?/ — sinx) (see Fig.l). One can see immediately that hiB) is a bounded 
function of two variables and hiA) is an unbounded. This means that quite different choice 
of function / is needed for these two cases in order to minimize the invariants. Inffuence of 
initial/boundary conditions is now also very clear - for instance, best approximation of hiB) 
can be obtained in the narrow belts of the lines parallel to one of the coordinate axis while for 
hiA) these directions are in no way special. 

To construct a sample of such an approximate factorization for the operator ((7j) we just 
suppose intuitively that auxiliary operator B is "good" if its coefficients differ from the coef- 
ficients of ((7j) not much, and its invariant is small. Our MATHEMATICA implementation of 
the BK-factorization includes simple graphic functions to display the differences between all 
the parameters of the initial and auxiliary operators. A choice of the function f{x, y) = sin ^ 

gives an auxiliary operator B of the form 

B = dxx — dyy + sin y sin — dx + cos x sin — dy -\ — (sin^ y — cos^ x) sin — . (8) 

xy xy 2 xy 

It is demonstrated at the Fig. 2 that for 10 < x,y < 100 operator B gives good enough 
approximation and correspondingly approximate factorization of the initial operator B has form 



B 



-( — cos X sin h sin — sin y^ + d^ + dy -( cos x sin h sin — sin y^ + d^ — 

2 V xy xy / J L2 V xy xy / 



xy xy 
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10-10 



10-10 



Figure 1: Invariant hi^) = \{,y^ — x"^) (left) and invariant h^B) = |(cos?/ — sinx) (right), in 
the domain —10 < x,y < 10 




100 100 



Figure 2: Upper panel: h^B) (left) and cio — aio (right); lower panel: aoi — aoi (left) and 
Ooo ~ (^00 (I'ight); in the domain 10 < x, y < 100 
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1 



Figure 3: Invariant hiB) in the domain 0.0001 < x,y <1 

with |/2(-B)| ~ 5 • 10^^. On the other hand, in the domain 0.001 < x,y < 1 quahtatively 
different approximation is needed while in this domain |/2(-B)| ~ 10^ (see Fig. 3). 

6 Brief Discussion 

We presented here an invariant formulation of BK-factorization which allows to factorize simul- 
taneously classes operators equivalent to the initial one under gauge transformations. We also 
showed the possibility to use the same procedure for the construction of the approximate fac- 
torization of LPDE in the case when corresponding LPDO is not factorizable. Obviously, if 
we get enough approximate factorizations of the given LPDE with different solvable first-order 
factors we can write out explicitly general solution of the initial LPDE. Otherwise, one gets a 
chain of the linear first-order equations 

AiO,n1pO = 0, Ao,n-l^l = V'O, •••• 

to be solved numerically which is a great numerical simplification, of course, specially for higher 
order LPDEs. On the other hand, while performing numerical simulations, one has to take into 
account a lot of other factors, first of all, initial and boundary conditions. It would be a 
nontrivial task to include them into the exact formulae given by BK-factorization. In order 
to estimate usefulness of this approach from numerical point of view we still have to answer 
all the questions concerning computation time, stability, computation error, etc. For instance, 
coming back to the example of approximate factorization given in the previous section, one 
have to estimate what is numerically more reasonable for a given set on initial and boundary 
conditions - to solve numerically the system of equations 

[i(cos a; sin ^ + sin ^ sin y) + - dy] o ipQ = 
[|(- cos xsmj^ + sin ^ sin y) + + dy] o = ipQ 

or one equation B oip = Q. Some answers can be given by the method presented in f5| where a 
symbolic approach is used to generate automatically finite difference schemes for LPDEs and to 
check their von Neumann stability. Some preliminary steps to be taken in this direction might be 
following: (1) to take a non-factorizable but solvable operator, for instance, Ai = dri.dy + xdx + 2, 
then LPDE Ai{ifj) = 0, has general solution 

i/j = -a^(^X(x)e-^'^ + f e'^y'^y^Y{y')dy''^ 
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with two arbitrary functions X{x) and Y{y); (2) to construct its approximate factorization 
Ai = Li o L2; (3) to get computational schemes using jSj - for Ai and Ai, (4) compute both 
numerically; (5) to compare results for Ai and Ai with the general solution for some classes of 
initial data and for a fixed choice of computational scheme. 
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